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ABSTRACT 

We analyze the possibilities of detection of hypothetical exoplanets in coorbital motion from 
synthetic radial velocity (RV) signals, taking into account different types of stable planar 
configurations, orbital eccentricities and mass ratios. For each nominal solution corresponding 
to small-amplitude oscillations around the periodic solution, we generate a series of synthetic 
RV curves mimicking the stellar motion around the barycenter of the system. We then fit the 
data sets obtained assuming three possible different orbital architectures: (a) two planets in 
coorbital motion, (b) two planets in a 2/1 mean-motion resonance (MMR), and (c) a single 
planet. We compare the resulting residuals and the estimated orbital parameters. 

For synthetic data sets covering only a few orbital periods, we find that the discrete radial 
velocity signal generated by a coorbital configuration could be easily confused with other 
configurations/systems, and in many cases the best orbital fit corresponds to either a single 
planet or two bodies in a 2/1 resonance. However, most of the incorrect identifications are 
associated to dynamically unstable solutions. 

We also compare the orbital parameters obtained with two different fitting strategies: 
a simultaneous fit of two planets and a nested multi-Keplerian model. We find that, even 
for data sets covering over ten orbital periods, the nested models can yield incorrect orbital 
configurations (sometimes close to fictitious mean-motion resonances) that are nevertheless 
dynamically stable and with orbital eccentricities lower than the correct nominal solutions. 

Finally, we discuss plausible mechanisms for the formation of coorbital configurations, 
by the interaction between two giant planets and an inner cavity in the gas disk. For equal 
mass planets, both Lagrangian and anti-Lagrangian configurations can be obtained from same 
initial condition depending on final time of integration. 

Key words: celestial mechanics, techniques: radial velocities, planets and satellites: forma- 
tion. 



1 INTRODUCTION 

The possible existence of exoplanets in coorbital motion has fas- 
cinated planetary scientists for several years. Since the diversity of 
exoplanetary configurations continue to surprise us, even more than 
15 years after the discovery of Peg5 lb, and it seems almost natural 
to expect Trojan planets to exist somewhere and the announcement 
of their discovery to be only a matter of time. Probably the first de- 
tailed analysis of hypothetical coorbital planets is due to Laughlin 
and Chambers (2002). They studied three types of coorbital con- 
figurations: tadpole orbits (around L4 and L5 equilateral points), 
horse-shoe configurations and "eccentric resonances". They pro- 
posed to distinguish between coorbital and single planet fits from 
RV data by observing residuals from long-term observations (more 
than ten orbital periods), because the coorbital configurations have 
large mutual interactions due to resonant motion. 

For systems with more than one planet, it is well known that 



the existence of resonant motion may be possible evidence of a 
past large-scale planetary migration due to interactions with the 
gaseous disk. Although their importance is unquestionable, it is still 
intriguing why some commensurabilities are very populated (e.g., 
2/1 MMR) and others that are currently empty (particularly the 1/1 
MMR). 

Lagrange (1873) discovered stable solutions for three massive 
bodies such that at all times their relative positions are located in 
the vertices of an equilateral triangle of variable size (L4 and L5 
solutions). Linear stability analyses traditionally focused on the re- 
stricted three-body problem (where one of the masses vanishes, e.g 
Morais 2001, Namouni et al. 1999 and references therein). Re- 
cently, Nauenberg (2002) numerically investigated the dynamical 
stability of general three body problem as a function of the eccen- 
tricity of the orbits and the Rouths' mass parameter. The resonant 
solutions were found stable when m ^j~ m2 < 0.038 and up to ec- 
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centricities ~ 0.6. The existence of such periodic orbits in the gen- 
eral three-body problem also was derived for example in Siegel 
& Moser (1971) and Laughlin & Chambers (2002). As mentioned 
in Laughlin & Chambers (2002), horse-shoe orbits are stable for 
planetary masses up to ~ 0Amj up , and despite first impression, 
Keplerian fit are good enough to reveal their radial velocity signa- 
ture with only four orbital periods. We note also that those type of 
orbits named "eccentric resonances", are simply an effect of an- 
gular momentum conservation in a coplanar three-body problem, 
that predicts when the eccentricity of one planet is maximum the 
other eccentricity reach a minimum and they have high amplitude 
of oscillation of resonant angles. 

More recently Hadjidemetriou et al. (2009) studied the prop- 
erties of motion close to a periodic orbit by computing the Poincare 
map on a surface of section. They constructed the symmetric fami- 
lies of stable and unstable motion describing one previous unknown 
symmetric configuration in the general problem, the Quasi-Satellite 
(QS) solution. 

In a previous work (Giuppone et al. 2010), we studied the sta- 
bility regions and families of periodic orbits of two planets in the 
vicinity of a 1/1 MMR, using numerical integrations and devel- 
oping a semi-analytical method. We considered different ratios of 
planetary masses and orbital eccentricities, although we assumed 
that both planets share the same orbital plane (coplanar motion). 
As result we identified two separate regions of stability, symmetric 
and asymmetric types of motion easily described with the behavior 
of resonant angles (a, Avj) = (A2— Ai, V02 — cci), summarizing] 

• Quasi-Satellite region: correspond to oscillations around a 
fixed point located at (a, Avj) = (0, 180°) independent of mass 
ratio. Although not present for quasi-circular trajectories, they fill a 
considerable portion of the phase space in the case of moderate to 
high eccentricities. 

• Lagrangian region: Two distinct types of asymmetric peri- 
odic orbits exists in which both a and Avj oscillate around values 
different from or 180°. The first is the classical equilateral La- 
grange solution associated to local maxima of the averaged Hamil- 
tonian function. Independently of the mass ratio mijm\ and their 
eccentricities, these solutions are always located at (a, Avj) = 
(±60°, ±60°). However, the size of the stable domain decreases 
rapidly for increasing eccentricities, being practically negligible for 
e, > 0.7. 

The second type of asymmetric solutions correspond to local 
minimum of the averaged Hamiltonian function. They were named 
Anti-Lagrangian points (AL4, and AL$) and, for low eccentrici- 
ties, are located at (a, Azu) — (±60°, =Fl20°). Each is connected 
to the classical Li solution through the cr-family of periodic orbits 
in the averaged system (a family of solutions where the angle a 
has zero amplitude of oscillation, e.g. Michtchenko et al. 2008ab). 
Contrary to the classical equilateral Lagrange solution, their loca- 
tion in the plane (a, Avj) varies with the planetary mass ratio and 
eccentricities. Although their stability domain also shrinks for in- 
creasing values of a, they do so at a slower rate than the classical 
Lagrangian solutions, and are still appreciable for eccentricities are 
high as ~ 0.7. An empirical relation exist for eccentricities below 
0.6 (e.g. Giuppone et al. 2010, Hadjidemetriou & Voyatzis 201 1) 



1 We recommend see the Figs. 2 and 5 from Giuppone et al. (2010) in 
order to easily identify these kind of coorbital motions in same plane of 
astrocentric parameters. We center this work on L4 and AL4 configurations 
since, due to symmetry consideration, they give same results as L5 and 
AL5 respectively. 
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In this paper we address the question of extrasolar Trojans, 
both from its detectability and possible origin. In a recent work, 
Anglada-Escude et al. (2010) showed that two planets in a 2/1 
MMR resonance could mimic a single planet in a more eccentric 
orbit. We wonder whether the opposite could also occur, and if two 
planets in a coorbital configuration could give similar RV signals 
as other dynamical systems including the 2/1 resonance. In Section 
2 we review the radial velocity signal generated by such systems 
and how they can be easily mistaken by the signal generated by a 
single planet. In the following section we select representative con- 
figurations and generate synthetics data-sets for several mass ratio. 
Sections 4 analyzes the dynamical stability of the different possi- 
ble configurations, while the dispersion of the best-fit parameters is 
discussed in Section 5. A comparison between two different strate- 
gies for planetary detection (two-planet fit vs. nested Keplerian fits) 
is presented in Section 6. Finally, a possible formation mechanism 
of extrasolar coorbitals is discussed in Section 7, where we present 
a series of numerical simulations (both hydro and N-body) of two 
planet systems immersed in a gaseous disk with an inner gap. Con- 
clusions close the paper in Section 8. 



2 KEPLERIAN RADIAL VELOCITY EQUATIONS 

First we need to understand some aspects involved in the RV sig- 
nal produced by two planets in coorbital motions, assuming unper- 
turbed Keplerian motion. Consider two planets with masses mi and 
m.2 in coplanar orbits around a star with mass mo = Mq . Let a; 
denote the semimajor axes, the eccentricities, Xi the mean lon- 
gitudes, vji the longitudes of pericenter, Mi the mean anomaly and 
fi the true anomalies. All orbital elements are assumed astrocentric 
and osculating. 

The suitable angular variables for coorbital motion are defined 
by g = A2 — Ai and Avj = VJ2 — vj\ . Disregarding mutual interac- 
tions between the planets, the radial velocity of the star is the sum 
of the individual Keplerian contributions (e.g. Beauge et al. 2007). 
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Ii being the orbital inclination with respect to the sky. Even though 
the radial velocity of a star perturbed by the two planets near the 
periodic orbit (i.e QS, Li and ALi) is not given by a single periodic 
signal, we may ask under what circumstances this signal could be 
mimic by a single planet orbiting the star or even two planets in 
other configurations. 

To address this question we can use the same approach than 
Anglada-Escude et al. (2010) rearranging the Keplerian contribu- 
tions from each planet. Using the expansions to first order in ec- 
centricities for the true anomaly / (see e.g. Murray & Dermott, 
1999) 

sin(/i) = sin(Mi) + e l sin(2Mi) 

cos(/ ! ) = cos(Mi) + ej cos(2Mj) - a, (3) 
and expanding 

cos (fi + vji) — cos (vji + Mi) + ei cos (2Mi + vji) — ei cos vji 
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Figure 1. Schematic configuration for coorbitals and corresponding single 
planet that mimic their RV signal. The position of single planet is deter- 
mined by equations (9) and (10} 



recalling that Ai = zui + Mi and substituting {4]l into ((2} we finally 
obtain 



V r 



K\ I cos Ai + ei cos (2Ai — vji 



+ K2 cos A2 + e2 cos (2A2 — ZU2] 



(5) 



Now, we suppose that the previous signal can be mimicked by 
a single planet. Imagine that its mean longitude A is located some- 
where between the position of coorbitals. According to Figure[T]we 
can rewrite the mean longitudes and the longitudes of pericenters 
as 



Ai = A — ai 

A2 = A + «2 , where 

zui = zu — j3± 

zu2 = zu + P2 where 



(6) 



ft + ft = Aro 



By definition of periodic orbit a = 0, Azb — giving con- 
stant values for ai, 0:2, ft and $2- We can then substitute eqs. {6} 
into eq. {5}, and regroup terms to rewrite 

V r = (Ki cos ai + K2 cos 02) cos A + e cos (2A — zu) (7) 
=K x 7 



K\e\ cos (2ai — ft) + K2C2 cos (2a 2 — ft 



(8) 



Kl COS Ql + 7\"2 COS «2 

where the conditions 

-fTisinai = Jf 2 sina2 (9) 
isTieisin(2ai-/3i) = K 2 e 2 sin (2a 2 - P2), (10) 

eliminate the additional periodic term and also define the position 
of fictitious single planet, giving 

K2 sin a 



tancti = 



tan(2ai-/?i) = 



#1 (1 + i|a C os a) 
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Figure 2. Parameters for a single planet to mimicking the radial velocity 
signal of two coorbitals with mass ratio 1712/ m\. The color code is used to 
distinguish configurations: QS in black, L4 in red and AL4 in green. Top: 
Values of Qi are shown in continuous lines while /3i is plotted in dashed 
lines. Middle: Associated value of K in units of K\. Bottom: Equivalent 
eccentricity of single planet for each configuration in units of the eccentric- 
ity of the coorbital planet with signal K\ . 



Qi = Q2 and fix — 02, from which the fictitious single planet 



would have A 



Ai+A 2 



and zu 



In this case ei = e2 



2 2 
and the RV signal arising from the single planet would be 

K = 2^i cos ^ 



ei- 



cos (a-^f) 



(12) 



(13) 



It can easily be seen that when both planets have the same mass 



In the opposite limit, when m,2 3> mi, both a.2^2 — > 
and the single planet is located near the position of the massive 
planet. Recall that for coorbital configuration with eccentricities up 
to et < 0.6, the stable zero-amplitude solutions satisfy eq. ([TJ, 
except for the L4 configuration that always lies in the line segment 
ei = ei. 

For other mass ratios, Figure|2]shows the equilibrium values of 
ati, Pi and the corresponding signal K amplitude and eccentricity 
e (the eccentricity is shown in units of ei). Note that the value of e 
only varies when the single planet mimics coorbitals in AL4. It is 
interesting remark that QS configuration could easily mimicked by 
single planet of circular orbit. 



3 ORBITAL FITS AND SYNTHETIC RADIAL 
VELOCITY 

First we review the differences in radial velocities signals using the 
Keplerian approximation l(2j and N-body integrator, in order to dis- 
tinguish how important are the mutual interactions and the feasibil- 
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Table 1. Arbitrary conditions near the stable periodic solutions in the 
(cr, Ara) plane, calculated with semi-analytical method. All conditions 
have a\ = ai- 



m2/ml = 1 



QS configuration 
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ity of being detected. The top frame of Figure [5] shows three syn- 
thetic curves covering four orbital periods for equally mass plan- 
ets (rrii — mj U p) in QS, L4 and AL4, configurations. The bottom 
frame of the same figure shows the error that would be obtained if 
the RV signal was calculated assuming constant Keplerian orbits. 
Even after four orbital periods, the error remains below ~ 4 m/s. 
As expected, the mutual perturbations would be even smaller if the 
individual mass were reduced. 

To test the mimicking effect described in the previous section, 
we selected six nominal solutions near exact fixed points chosen 
from Giuppone et al. (2010), three with mi = 7712 and three other 
with mi = 7712/3. In all cases we fixed 7712 = lmj up . Initial con- 
ditions are summarized in TableQ] In order to test the validity of our 
expressions l[5) and {7]l-l[8]l, constructed from a first-order expan- 
sion in the eccentricities, here we deliberately chose high eccentric 
orbits (ei = 0.3). If our predictions prove correct in these cases, 
we can be assured that they will be valid for lower eccentricities as 
well. 

For each nominal configuration we generated a synthetic RV 
curve describing the stellar motion around the barycenter of the sys- 
tem. The curve was then represented with a discrete sampling of TV 
observation times U distributed randomly, according to an homoge- 
neous distribution (thus avoiding aliases in data from daily/seasonal 
observations, e.g. Dawson & Fabrycky, 2010). In each data point 
the nominal RV value V r {ti) was displaced to a new value V ri — 
V r (ti) + A/"(0, e) following a Gaussian distribution with constant 
variance e 2 . The resulting synthetic data set was the used as input 
for our orbital fitting code PISA (Pikaia genetic algorithm + simu- 
lated annealing) (e.g. Beauge et al. 2008). Since our data sets cover 
only four orbital periods, the orbital fit was obtained assuming non- 
interacting Keplerian orbits. 

1000 different data sets were generated for each orbital con- 
figuration (QS, 1/4, AL4) and each mass ratio 7712/7711. Each syn- 
thetic data set consisted of N — 100 data points covering a total of 
four orbital periods. The results shown here were obtained adopting 
e = 3 m/s, a value similar to the present day errors in detections 
programs (including stellar jitter). The same analysis was repeated 
for other values of e, showing similar results although the disper- 
sion around the mean values decreased/increased as function of e 
(see e.g. Giuppone et al. 2009). 

As noticed by Laughlin & Chambers (2002) and Gozdziewski 
& Konacki( 2006), two planets in a coorbital configuration produce 
only one peak in a Fourier spectrum, meaning that it could easily 
be confused with a single planet (obviously more massive than the 
individual components). On the other hand, Anglada-Escude et al. 
(2011) found that two planets in circular orbits near a 2/1 MMR 
could also be falsely detected as a single planet in an eccentric orbit. 
So, our question here is the following: is it possible for two planets 
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Figure 3. Top: Synthetic radial velocity curves generated with an N-body 
integrator for three different orbital configurations, specified in the top of 
the graph for the first conditions in Table [T] Bottom: Difference between 
radial velocities calculated from the N-body integration and those gener- 
ated with Keplerian model. Although the difference increases with time, it 
remain below 4 m/s even after four orbital periods. 

in different coorbital configurations (QS, L4 or AL4) to be falsely 
identified as a single planet or as two bodies in a 2/1 resonance? 

To test this idea, each of the synthetic data sets was fitted as- 
suming: (i) a single planet, (ii) two planets in coorbital configura- 
tion (without specifying the type) and (iii) two planets in the vicin- 
ity of a 2/1 MMR. The resulting values of wrms (weighted rms) for 
the 1000 data sets associated to each nominal solution are displayed 
as histograms in Figure|4] Each frame corresponds to a given nomi- 
nal solution and mass ratio, while the color code for the histograms 
corresponds to the fitted orbital configurations. For example, in the 
top left-hand plot we constructed 1000 synthetic data sets from 
a nominal quasi-satellite configuration of two planets with equal 
masses. The data sets where then fitted assuming the three different 
possible configurations, the resulting values of wrms plotted. The 
colored histograms show the distribution of the residuals of each 
assumed configuration. As we can see, the results corresponding to 
QS solutions show to smallest distribution of wrms (both in mean 
and dispersion), very similar to the adopted value of e. Recall that 
wrms ~ e is usually employed as evidence that not only the fit is 
satisfactory but also that the assumed planetary model is correct. 

Although for nominal QS solutions and mi = 7712 the coor- 
bital configuration gives the best fits, this is not always the case. 
For systems with a more massive outer body (7712 /mi = 3), and 
for any nominal coorbital motion, we find that all three assumed or- 
bital solutions given practically the same error distribution (right- 
hand side plots). In other words, it seems that for 7712 > mi the 
difference in short-term RV signals of one planet or two bodies in 
either a 1/1 or 2/1 MMR are virtually indistinguishable. For equal- 
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Figure 4. Histograms of residuals adopting different orbital configurations 
to fit a series of synthetic RV signals. The nominal solution is shown in 
the top left-hand side of each frame. Each colored histogram correspond to 
a different assumed orbital configuration: two planets in coorbital motion 
(black), a single planet (red) and two planets in a 2/1 MMR (green). 



Table 2. Solutions from synthetic data sets generated from configuration 
AZ/4 in Table 1 . 
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Figure 5. Comparison of synthetic curves (top frame) and residuals (bottom 
frame) from a nominal AL4 solution and fitting the synthetic data set with 
three planetary models: coorbitals (black), single planet (red) and 2/1 MR 
(green). The best fit values of the masses and orbital elements are shown in 
Tabled 



fit mi m,2 a± 0,2 e\ e2 wrms 

model nijup rn-jup AU AU m/s 

AL4 0.37 0.9 1 1 0.31 0.127 3.03 

1 Planet 1.12 - 1 0.107 - 3.26 

2/1 MMR 0.18 1.14 0.63 1 0.18 3.05 



mass planets, the picture is similar, specially if the nominal con- 
figuration corresponds to two planets in L4. However, two bodies 
in other coorbital configurations (i.e. QS or AL4) appear easier to 
identify. 

Figure [5] shows an example. The top frame shows a synthetic 
data set constructed from a nominal solution of two planets in 
a stable AL4 configuration. The three different fits are superim- 
posed in different colors and show practically no difference. The 
corresponding distribution of residuals are shown in the bottom 
frame, also displaying no significant difference. The best fit plan- 
etary masses and orbital elements of each orbital configuration are 
shown in Table 2. 



4 STABILITY ANALYSIS OF THE FITS 

Although a given coorbital pair may be falsely identified as another 
orbital configuration, we still do not know whether all possible al- 
ternatives are equally dynamically stable. To analyze this point, we 
numerical integrated each of the best fits for a time span covering 
10 4 orbital periods, and calculated the amplitude of oscillation of 



the resonant angle a as well as the difference in longitudes of peri- 
center Azu = ti72 — zo\ . In the case of coorbital motion the resonant 
critical angle is basically the synodic angle a = A2 — Ai, while for 
the 2/1 MMR we have chosen the so-called principal resonant angle 
<j = 2A2 — Ai — W\, 

The Table [3] presents a statistical analysis of the dynamical 
evolution of all the best fits. We denote the dynamical behavior as 
"resonant" if the system is stable during the full integration span 
and both a and Avj show small-to-moderate amplitudes of oscil- 
lation (< 50°) around the fixed point. We call librators those con- 
figurations where the resonant angle librates but the secular angle 
circulates. "Non-resonant" configurations are those for which both 
angles circulate. Finally, "unstable" configurations are those that 
led to collisions or escapes of one of the bodies within the integra- 
tion span. 

The left-hand side of Table [3] describes the dynamics of the 
coorbital fits. For both mass ratios, the QS solutions appear the 
most robust. For equal mass planets, all the fits correctly yield sta- 
ble resonant quasi-satellite orbits, fully consistent with the nominal 
solution. For 1712/ mi = 3, however, 13% of the synthetic data 
sets led to unstable coorbital solutions. The reliability of the fits 
decreases when the adopted nominal solution is L4 and even worse 
when it is AL4. In the latter case, more than a third of the data sets 
led to unstable solutions. 

The right-hand side of the Table 3 now shows results for those 
fits that incorrectly identified the signal as that generated by two 
planets in the vicinity of the 2/1 MMR. For equal mass planets, the 
best fits are highly unstable for both QS and L4 nominal configura- 
tions, while the AL4 data give a 63 % of stable configurations but 
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with both resonant angles circulating. When we increase the mass 
ratio, the compatible 2/1 MMR solutions are sightly more stable 
although again correspond mainly to non-resonant configurations. 
Again, the QS and L4 configurations are the most unfavorable for 
being labeled as 2/1 MMR solutions, although for asymmetric equi- 
lateral Lagrange solutions, a small fraction (~ 8%) of 2/1 MMR 
solutions show stable librators. 



5 DISPERSION OF THE BEST-FIT PARAMETERS 

Although the best fits assuming two planets in a 2/1 MMR have 
proved to be largely unstable, and thus difficult to confuse with 
coorbital planets, we have no clear indication of the cause of the 
instability. In this section we will analyze the planetary masses and 
eccentricities obtained form the diverse best fits, and compare them 
with the nominal solutions. 

Results are presented in Figures [6l8l for equal mass planets 
(mi = m2 = injup). Each corresponds to a different nominal 
coorbital solution (QS, L4 and AL4, respectively), and is divided 
into four frames. The top graphs show histograms with the dis- 
tribution of the eccentricities (left) and deduced planetary masses 
(right). Masses are in units of mj llp . Different colors and line 
types correspond to different configuration or bodies (see caption 
for details). The two bottom plots show the dispersion of the ec- 
centricities (left) and relation between the mass and eccentricities 
(right). Once again, different colors are used for different configu- 
rations/bodies. 

For a nominal QS configuration (Figure|6j, the best-fit param- 
eters assuming a coorbital solution show a considerable dispersion: 
approximately 0.04 in eccentricities and 0.1mj up in masses. How- 
ever, there appears to be no appreciable systematic error and both 
distributions seem symmetric with respect to the nominal values. 
The single-planet fits show a very steep distribution around e ~ 
and m ~ 2. To understand this result, we can analyze the expected 
RV signal. For a QS configuration we have a = and Azu = 180° . 
Moreover, from equation {T} we have that equal mass planets imply 
ei = e2, from which we expect Ki — K2. Consequently, equation 
leads to 



V r = 2Ki cos A. 



(14) 



This tells us that two planets in QS orbit produce the same signal as 
one planet with semi-amplitude K — 2K\ and eccentricity e = 0. 
The dispersion of the best fits show an excellent agreement with this 
prediction, yielding a planet with m = 7711+7712 in a quasi-circular 
orbit. Surprisingly, the dispersion of the different fits is smaller than 
obtained assuming (correct) QS orbits. 

The distribution of the 2/1 MMR fits, shown in the histograms 
with dashed lines, show a strong bimodal shape. Approximately 
half of the best fits give masses in the vicinity of (7711,7712) ~ 
(2,0.5) and eccentricities around (ei,e2) ~ (0,0.3). The other 
half of the solutions clutter around (7711,7712) ~ (2,0.1) and 
(ei,ea) ~ (0.1, 0.8). This indicates the existence of two local min- 
ima in the residual function with similar values of the wrms and 
similar extension in the parameter space. Small differences in the 
synthetic data sets would then highlight one or the other, thereby 
yielding either two planets in low-to-moderate eccentricities or a 
solution in which one of them is almost parabolic. However, as 
shown in Table[3] practically all solutions are unstable. 

Results for a nominal L4 configuration are shown in Figure 
[7] Contrary to the previous case, the coorbital fits now show a very 
narrow dispersion around the correct solution, implying that the RV 




Figure 6. Dispersion of best-fit solutions from synthetic data sets for two 
planets with mi = mi = 77ij up in a QS solution. Top: Distribution of the 
best-fit eccentricities (left) and planetary masses (right). The single-planet 
fit is shown in a continuous green line. For the two-planet fits, the results 
for planet 1 are shown in black, while those corresponding to planet 2 are 
shown in red. In these cases, continuous lines are used for the coorbital 
configuration, while dashed lines correspond to a 2/1 MMR fit. Bottom 
left: Scatter plot showing the dispersion of best-fit eccentricities. Results 
for the coorbital fit are shown in black, while those corresponding to a 2/1 
MM fit are shown in gray. Bottom right: Scatter plot with the dispersion of 
eccentricity as function of planetary mass. Black and red dots correspond 
to both planets of a coorbital solution, orange and gray correspond to the 
planets of a 2/1 MMR and green points are the single planet solutions. Gray 
circles were drawn to identify nominal configurations. 
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Figure 7. Same as previous figure, but for two planets with mi = 7712 in 
an L4 configuration. 
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Table 3. Percentage of different dynamical outcomes of the best fits for different mass ratios mi/mi. The value of mi was fixed at lmj up . See text for a 
detailed description. 
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signal of two planets in an equilateral configuration is more robust 
than two bodies in QS. Applying equation to this case yields 
Qi = Q2 = 30°, and the same RV signal could also be associated 
to a single planet with 



Vr 



V3Ki\ cos A + e cos (2A — zu) 



(15) 



where e = e\ = e-2 is the value of the eccentricity of the original 
bodies. The distribution of the best single-planet fits follow this 
prediction and show a very sharp peak at (m, e) = (~ 1.7, 0.3). 

The 2/1 MMR best fits always give two planets with a large 
mass ratio. The largest body has mass and orbital elements similar 
to the single planet fit, while its companion has very small mass 
and very high eccentricity (K — > a and e — > 1). Thus, this second 
planet does little more than attempt to resolve the residuals of the 
single planet solution. It is therefore no surprise that practically all 
resonant configurations are highly unstable. 

Finally, Figure [8] shows the results obtained from a nom- 
inal AL4 configuration. Initial conditions show a small ampli- 
tude oscillation around the equilibrium solution (a, Azu) — 
(95.4°, 257.46°) which, substituting into equation 10, leads to an 
equivalent single planet solution with 



Vr 



l.3AKi [cos A + 1.24 ei cos (2A - to)] , 



(16) 



which results in the strong peak in the mass and eccentricity his- 
tograms. However, all one-planet solutions have much larger resid- 
uals than the coorbital solution (see Figure[4]l. Once again, the so- 
lutions assuming a 2/1 MMR generate a bimodal distribution, with 
part of the solutions (~ 33%) with high eccentricities and conse- 
quently unstable. However, contrary to the previous cases, almost 
two thirds of the fits lead to configurations that are dynamically sta- 
ble, even if the resonant angles are in circulation. These correspond 
to mi ~ 0.4 and ei < 0.2. 

The same procedure was also followed for synthetic data sets 
with mj/tni = 3 and ei = ei/3. Results showed very similar 
traits as those presented for equal-mass planets, except for a larger 
dispersion in the parameters. This larger dispersion diminishes the 
proportion of stable coorbital fits (see Table [3]( but also allows for 
a significant number of stable fits assuming two planets in a 2/1 
MMR. In the case of a nominal AL4 configuration, almost one third 
of the synthetic data sets gave stable 2/1 configurations, although 
only a fraction of these corresponded to libration. 
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Figure 8. Histograms of fitted synthetic data sets for two equally mass 
planet in AL4 configuration. 



6 NESTED VS. TWO-PLANET FITS 

Apart from possible misidentification of coorbital planets as other 
configurations, the 1 : 1 MMR also appears sensitive to the strategy 
adopted for the fitting algorithm. There are two commonly used 
procedures to fit a two-planet solution into a given RV data set. 
One possibility, the simultaneous two-planet fit, assumes the exis- 
tence of two masses from the beginning, and attempts to fit the data 
set with a model with two periodic signals. The second approach, 
sometimes referred to as the nested model first attempts to fit one 
planet into the data. If the residuals are too large or if its Fourier 
spectrum shows a significant periodicity, then a second planet is 
fitted into the reduced data. This typically is adopted if the largest 
amplitude of the spectrum is larger than a given false alarm proba- 
bility (FAP). 

In a perfect world, both procedures should give the same re- 
sults, or at least very similar to each other. This appears to be the 
case for non-resonant planets or even for planets in the 2/1 MMR 
(see the analysis done by Beauge et al 2008 for the HD82943 sys- 
tem). However, as we will show below, they can lead to completely 
different solutions in the case of coorbital planets. 

Using the same type of nominal configurations shown before, 
we generated new synthetic RV data sets. We incorporated two dif- 
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Table 4. Best multi-Keplerian orbital fits for the data generated by two equal 
mass planets in a QS coorbital configuration. The time of passage through 
the pericenter r is given with respect to the time of the first data point. All 
orbital elements as astrocentric and osculating. 





Simultaneous 2-Planet Fit 


Nested 2-Planet Fit 


Parameter 


Planet 1 


Planet 2 


Planet 1 Planet 2 


K [m/s] 


28.59 


31.00 


54.19 5.26 


P[d] 


396.9 


402.2 


399.62 132.25 


e 


0.3421 


0.2696 


0.0194 0.1835 


m [mj up ] 


1.8757 


1.8469 


1.96 0.13 


a [AU] 


0.7461 


1.1779 


1.06 0.508 


V [m/s] 


-0.051 




0.05 


wrms [m/s] 


5.004 




5.62 




1.027 




1.154 



ferences: first, the orbital period of the planets was raised to 400 
days. Second, the data covered 12 years (10 complete periods) with 
N — 200 randomly distributed observations assuming e=5 m/s. We 
repeated the same test fits as in the previous section with similar re- 
sults. 

We now use the same data to compare both fitting strategies. 
Results are shown in Figure [9] All the periodograms were calcu- 
lated using DCDFT (Ferraz-Mello 1981) which allows a treatment 
of unequally spaced data. In the top and middle frames the power 
spectra are normalized so that the total area is unity. The top graph 
shows the periodogram of the original data, where the main 400 day 
signal is clearly visible. The middle frame shows the power spec- 
tra of the residuals after a single planet fit. The dashed horizontal 
line corresponds to a false Alarm Probability of FAP = 10 -4 , and 
was estimated using the Quast algorithm (Ferraz-Mello & Quast 
1987). We note the existence of peak corresponding to a period 
near P = i P% which appears statistically relevant. After perform- 
ing a second planet fit the wrms was reduced from 6.6 m/s to 5.6 
m/s. 

Finally, the lower frame shows periodograms of the residuals 
after both two-planet fits. Black curves show the results employ- 
ing a nested model, while the red curve correspond to a simultane- 
ous two-planet model. The power spectra are quite similar and no 
additional periodicity is observed. The resulting masses and plan- 
etary parameters of both models are given in Table [4] Although 
both results are very different, not only are the wrms very similar, 
but the incorrect system derived from the nested model is also dy- 
namically stable. In fact, the orbital eccentricities derived from the 
nested model are actually significantly lower than the nominal val- 
ues. Finally, we found no appreciable difference by using N-body 
fits for this data set. 

The top frame of Figure [TOl shows, in blue, the data synthetic 
points of the fictitious QS configuration. The black curve presents 
the expected RV signal from a nested two-planet fit, while the red 
curve shows corresponds to a simultaneous fit of both planets. The 
bottom plot gives the resulting o — c showing practically no differ- 
ence between both fits. 



7 ORIGIN OF EXOTROJANS 

In the first part of this paper we have analyzed whether a coor- 
bital configuration could give RV signals similar to other dynami- 
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Figure 9. Top: Power spectra of a synthetic data of two coorbital planets 
with orbital period of 400 days. Middle: Power spectra of the residuals after 
a single-planet fit. Horizontal dashed lines marks the 99.99% confidence 
level. Bottom: Power spectra obtained from the residuals of two different 
two-planet fits. Red correspond to a simultaneous fit of two massive bodies, 
while the black curve is the results of a nested two-planet model. 



cal systems, which may explain why so far no coorbitals have been 
discovered. In this second part we focus on a possible formation 
mechanism. 

There is a vast literature trying to explain the formation of 
Trojan planets. Several different mechanisms have been proposed, 
including accretion from protoplanetary disc (Laughlin & Cham- 
bers 2002) pull-down capture into the 1/1 resonance, direct col- 
lisional emplacement, and in situ accretion (Chiang & Lithwick 
2005), or convergent migration of multiple protoplanets (Thommes 
2005, Cresswell & Nelson 2006). 

Kortenkamp (2005) studied the gravitational scattering of 
planetesimals by a protoplanet, revealing that a significant fraction 
of scattered planetesimals (between 5 to 20%) can become trapped 
as QS in heliocentric 1/1 coorbital resonance with the protoplanet. 
They included a solar nebula gas drag and considered planetesi- 
mals with diameters ranging from ~1 to ~1000 km. The initial 
protoplanet eccentricities were chosen from e = to 0.15 and pro- 
toplanet masses range from 300 Earth-masses (7719) down to 0.1 

A different possibility was analyzed by Cresswell & Nelson 
(2006), where the authors used a hydrodynamical code to simulate 
the evolution of systems of up to 10 planets with masses between 
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time [days] 



Figure 10. Top: Synthetic data set of observations (marked as blue dots) 
and two curves generated with solutions obtained with coorbital simultane- 
ous fit (red) and two nested single fits (black). Bottom: Residuals from both 
fits showing no significant differences. 



2m@ and 20m^ . After brief period of chaotic interaction charac- 
terized by scattering, orbital exchanges and collisions, some cases 
led to coorbital planets occupying either horse-shoe or tadpole or- 
bits. If the initial separation between bodies was taken A ~ 5i?mn 
such configurations occurred in 20% of the runs, while the prob- 
ability increased to 80% if the initial separation was reduced to 
A = 4i?Hiii- They also found that planets sometimes evolved from 
horse-shoe to tadpole configurations, although some coorbital con- 
figurations were shot-lived, ultimately leading to disruption and es- 
capes of one of its members. 

Beauge et al. (2007) studied the in-situ formation of 
terrestrial-like Trojan planets starting from a swarm of planetesi- 
mals initially located around the Lagrangian point L4 of a massive 
planet. They analyzed different masses, initial conditions and both 
gas-free and gas-rich scenarios. The swarm of proto-trojans were 
initially located within the L4 tadpole region with total masses 
ranging from 1 to 3me. The gas interaction was modeled using 
linear and non-linear gas drags in a N-body code, adjusting the 
parameters to reproduce preliminary tests done with the FARGO 
hydrodynamical code. Their main results showed that the accre- 
tional process within the coorbital region is not very efficient, and 
the mass of the final Trojan planet never seems to exceed 0.6 m^g. 

Hadjidemetriou & Voyatzis (2011) studied the evolution of a 
QS system under the additional effects of a gas drag. They found 
that a QS planetary system with initially large eccentricities can 
migrate along the family of periodic orbits and be finally trapped 
in a satellite-type orbit. Thus the authors provided a mechanism for 
the generation of satellite systems starting from a planetary config- 
uration. These results agreed with those obtained numerically by 
Kortenkamp (2005). 

Morbidelli et al. (2008) complete this brief picture, analyz- 
ing the interactions between several planetary cores and a "density 
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Figure 11. Surface density of the gas disk in the vicinity of a density jump 
or inner cavity. E D and E; are the values outside and inside the jump, also 
characterized by ahalf-width A. The black curve shows the azimuthal aver- 
age obtained from a hydrodynamical simulation (FARGO) while the 
red curve is the analytical approximation using equation 117t . 

jump" in a gaseous disk. Depending on the planetary masses and 
initial conditions, they found cases of resonance trapping, scatter- 
ing and even the formation temporary binaries. 

In this paper we follow the same scenario of Morbidelli et al. 
(2008) to analyze whether coorbital systems may also be formed 
through the interaction of two planets with a density jump in pro- 
toplanetary disk. We thus consider two planets initially far from 
mean-motion commensurabilities that migrate inwards, ultimately 
becoming locked in resonance. We assume that the disk possesses 
an inner cavity which acts as a planetary trap, halting the migration 
of the massive bodies. For simplicity, we consider a constant sur- 
face density S outside the cavity, and a constant value Si inside. 
The density jump is set at a radial distance r = 1 (arbitrary units) 
from the star. Thus, we characterize the density jump in the disk by 
two main parameters: the density ratio F — E /Ei and the width 
A of the cavity edge (see Figure [Til. 

Our main simulations were performed using the FARGO-2D 
public hydro-code, whose basic algorithm may be found in Masset 
(2000). This code solves the Navier-Stokes equations for a Kep- 
lerian disk subject to the gravity of the central object and that of 
embedded planets. The disk is assumed isothermal and not self- 
gravitating. 

7.1 Preliminary N-body simulations 

Before undertaking the main hydro-simulations, we wish to ana- 
lyze what disk properties (S , F, A, etc.) and planetary parame- 
ters (masses, initial separation, etc.) are compatible with the for- 
mation of coorbital solutions. Since such a broad sweeping of the 
parameter space is extremely time-consuming for a hydro-code, we 
employed an N-body code for a series of preliminary runs. The in- 
teraction between the planets and the disk was approximated by a 
Type I migration following the semi-analytical model of Tanaka et 
al. (2002) and Tanaka and Ward (2004). In the case of a smooth 
density profile, explicit expressions for the forces acting on the 
planetary masses can be found in Ogihara et al. (2010). 

Since we will be working with planetary masses smaller than 
one Saturn-mass, Type I migration is expected to be a fair descrip- 
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tion of the tidal interactions with the disk. However, we stress that 
these N-body simulations are not expected to be accurate, but are 
used solely as guidelines for the full hydro-simulations introduced 
in the following sections. 

The inner cavity in the disk may be approximated with an hy- 
perbolic tangent, such that 

r- 1 



4 -i 



In E(r ) = C a tanh x + Cb with 
and where the coefficients are given by 



C a 



2 (l n E o 



•InEi) 



C b = -(In S +lnE i; ) 



(17) 



(18) 



Figure [TT] shows two representations of S(r) with a test cavity 
(edge at r — 1). The black curve was generated with FARGO 
(see Benitez-Llambay et al. 2011 for details), while the red curve 
was obtained applying equation J17| >. Both appear very similar, al- 
though the analytical function has a smoother trend near the edges 
of the cavity. 

Having an explicit functional form for the density profile, we 
can modify Tanaka's equations for the corotational torque to re- 
produce the contribution generated by a density jump in the disk. 
This was done following the calculations deduced in Masset et al. 
(2006a), and the resulting expressions were incorporated into our 
N-body code. With this tool we were then able to reproduce the dy- 
namical behavior of massive planets near the density jump without 
needing to artificially halt the planet at a given orbital distance from 
the star. 

Benitez-Llambay et al. (201 1) showed that the real population 
of close-in exoplanets is consistent with a disk inner cavity located 
at ro ~ 0.01 AU. However, since our N-body simulations will 
be used primarily as a stepping stone into hydro-simulations, it is 
preferable to scale the initial conditions to ro = 1. The Appendix 
shows the scale transformations necessary for the gas density and 
time units to ensure the same dynamics. For example, assuming a 
typical surface density for the disk of E ~ 5800 gr/cm 2 at r = 1 
AU, the corresponding value at any arbitrary spatial scaling param- 
eter ro gives E = 6.4 x 10 -4 MqTq 2 , where ro is given in AU. 

Figure [T2l shows the results of several simulations using our 
modified N-body code with a density jump centered at r = ro = 1 
and with a half-width of A = 2H/r. Throughout all our simu- 
lations, both N-body and hydro, we assumed a disk scale-height 
equal to H/r — 0.05 and a depth of the cavity equal to F = 
E /Ei = 10. The density outside the cavity will be set to E = 
TVEo with N a positive integer and E = 10" 3 M Q AU~ 2 . In 
other words, we will perform simulations with different surface 
densities given as multiples of an arbitrary base value Eo. 

All plots shows the time evolution of the ratio of mean- 
motions ni /712 between both planets. In the top plot we considered 
a fixed value of E = 15Eo, but varied the mass ratio of the bodies 
in the interval 7712/mi £ [10~ 2 , 10°]. In all cases, the mass of the 
largest body was equal to one Saturn mass. 

Since mi has the smallest initial semimajor axis, as well as the 
largest mass, its orbit suffers a faster orbital decay, and is the first 
planet to reach the cavity edge. In all the runs, this body was trapped 
near the center of the density jump (r ~ 1.05ro) with a small ec- 
centricity. Meanwhile, the outer approaches the cavity edge, en- 
countering successive and increasingly strong mean-motion reso- 
nances in its path. The orbital evolution is then dictated by two 
opposing forces: (i) the differential disk torque that subtracts angu- 
lar momentum and causes further orbital decay, and (ii) resonant 
perturbations with the inner mass that can generates conditions for 
a net increase in the angular momentum. 
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Figure 12. N-body simulations of two planets encountering an inner cav- 
ity in the disk located at r = ro = 1 (arbitrary units). Top: Using constant 
value for density (S = 15 X So and different pairs of planets with mass ra- 
tio from 10 -2 to 1. The coorbital configuration is achieved if the mass ratio 
is 0.7 < 1712/mi < 1-4, this range becoming wider for increasing values 
of E . Middle Frame. Two Saturn-mass planets in a disk with different E Q . 
For values larger than 12Eo final coorbital configurations appear with high 
probability. Bottom Frame. Two simulations with mi = m,2 = msat 
and E = 15So, but considering different magnitudes for the eccentricity 
damping. Meanwhile planet 1 was fixed at 1.5 AU, we change the position 
of planet 2 from 1.9 to 3.2 AU in order to represent better the results in the 
figure. 



The top plot of Figure Q2] shows that the system is stopped 
at different mean-motion resonances depending on the mass ra- 
tio. For small values of m^jmi the two planets are trapped in the 
3/1 MMR, while the 3/2 commensurability is favored for larger 
mass ratios. This is expected, since a larger outer mass requires 
stronger resonance perturbations to counteract the disk torque. For 
m2/mi > 0.7, however, it appears no commensurability is suf- 
ficiently strong, and both bodies reach the cavity edge and evolve 
towards a stable coorbital configuration. 

From these results it then appears that coorbital planets may be 
formed at a density jump in the disk if the mass ratio is sufficiently 
close to unity. We checked that the same result was obtained for 
other mass values, considering planets ranging from Neptune to 
Jupiter analogues. 

The middle plot of Figurell2lnow analyses the dependence on 
the value of the surface density E outside the cavity, while main- 
taining the same cavity depth F. Results are shown for two equal 
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Figure 13. Comparison between two hydro-simulation with the same ini- 
tial conditions (see text for details) but different grid resolutions. The high 
resolution simulation was constructed with 306 radial zones, while the low 
resolution run contained only 168 radial zones. Both lead to coorbital con- 
figurations with similar global dynamical features 




mass planets (mi = m.2 = ms a t) and four different multiples of 
the base density Eo. We note that smaller densities favor resonance 
trapping in non-coorbital configurations (e.g. 5/4 MMR), but for 
large values of E , typically above Ei ~ 10Eo, coorbital systems 
are the usual outcome. Of course there is some dependence on the 
initial conditions, mainly on the initial separation of the bodies, but 
these results are surprisingly robust and are representative of most 
of our runs. 

Finally, the bottom frame shows the dependence on the ec- 
centricity damping force Fdamp- Although both the circularization 
time scale r e and the orbital decay time r a have the same depen- 
dence with the planetary mass and disk density (e.g. Tanaka and 
Ward 2004), it is possible to fictitiously modify the eccentricity 
damping by changing the tangential force suffered by the bodies. 
Our main aim in this experiment is to see how the dynamics may 
be affected by different values of r e /r a . The plot shows two sim- 
ulations, one with the correct value of r e (black) and one in which 
the damping rate was reduced by a factor 10 (blue). Although a 
coorbital configuration was achieved in both cases, the reduced 
damping generates a much larger amplitude of oscillation of the ec- 
centricities which may ultimately lead to orbital instability. In fact, 
some of these simulations showed complex exchanges between dif- 
ferent types of coorbital solutions (L4, L5, etc). 

In conclusion, our preliminary N-body simulations appear to 
indicate that coorbital systems of massive planets can be obtained if 
the mass ratio is sufficiently close to unity and the surface density 
of the disk is sufficiently high. Although the final outcomes have 
a certain dependence on the initial conditions (like all resonance 
trapping phenomena), we repeated the runs for initial separations 
between 0.5 to 1.7ro, finding no significant changes in the results. 
Inspired by these results, our next step is to upgrade to hydrody- 
namical simulations. 



7.2 Hydro-numerical setup 

Our hydro-simulations were carried out using the FARGO code 
(Masset, 2000) with a central star of one solar mass. We consider 
two-dimensional, locally isothermal disks uniformly distributed 
over a polar grid. 

The disk inner cavity (IC) was generated using an ad hoc step 
in kinematic viscosity around r = ra = 1, using the same recipe 
as described in Masset et al. (2006a). For all our runs we adopted a 
cavity depth of F = E /Ej = 10, and a half-width of A = 0.4. 
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Figure 14. Snapshot of the final outcome of a high resolution (306 radial 
zones) hydro-simulation leading to two planets locked in a coorbital con- 
figuration. One of the planets (mi) is fixed in the x-axis. The numbers 1, 2, 
3 and 4 represent the location of the other planet (pi2) at different times of 
the run, with number 4 being the final configuration. 



The kinematic viscosity was considered uniform outside the cavity, 
and equal to v = 10 ro^o , where fio is the orbital frequency 
at radius rn. The disk aspect ratio was chosen as H/r = 0.05 and 
also constant over r. Boundary conditions were chosen to be non- 
reflecting for the inner edge of the disk. 

The disk was represented with a spatial resolution of 384 az- 
imuthal zones and 168 radial zones. The inner radius of disk is 
set to 0.42ro and the outer radius is set to 4.2ro. This resolution 
allows us to identify details as small as 0.04ro. Casoli & Masset 
(2009) showed that the coorbital co-rotation torque exerted by the 
disk over a planet occurs in a scale of x s ~ O.Olro for m ~ 2m®, 
while Masset et al. (2006b) showed that the horse-shoe zone width 
in a two-dimentional disk scale as x s = 1.16^/ q/h. Consequently, 
x 3 should increase with the planetary mass. For example if we are 
interested in the dynamical evolution of Saturn-type planets, we 
expect that our choice for the spatial resolution should be adequate. 

To test this idea, we performed several simulations with the 
same initial conditions but different spatial resolutions for the disk, 
and compared the dynamical evolution of the planets. An example 
is shown in Figure[T3"1for two Satum-like planets initially in circu- 
lar orbits far from the cavity edge and in non-resonant conditions. 
Both simulations show the same qualitative results, not only in the 
final outcome of the planets (both runs end in coorbital configura- 
tions) but also in the behavior of the eccentricities. The only notice- 
able change is a slight difference in the radial distance at which the 
planets are trapped. However, this is not significant and could be 
explained due to a decrease of the co-rotational torque from a low 
resolution description of the horse-shoe region. 

Following the predictions of our N-body runs, all our hydro- 
simulations were done with two equal-mass planets (mi = m^) 
in a high-density disk. Initial conditions were varied both with re- 
spect to the initial separation between the planets and the initial ec- 
centricities. We also considered different planetary masses between 
one Earth-mass and one Satum-mass. 
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Figure 15. Results of a low resolution (168 radial zones) hydro-simulation 
of eight planets with masses in the Earth-Saturn range. Top plot shows the 
semimajor axis as function of time, while the lower frame shows the eccen- 
tricities. 



Under a wide variety of initial conditions, we found that the 
orbital evolution of two Saturn-mass planets ultimately lead to the 
formation of stable coorbital configurations. An example is shown 
in Figure[l4] where both bodies are trapped in Li/ Xs-type motion 
with a libration of both the resonant angle a = A2 — Ai and the 
difference in longitudes of pericenter Azu = ZU2 — zji . 

To test the robustness and stability of the coorbital solutions to 
additional perturbations, we performed a new series of runs includ- 
ing additional six Earth-mass planets initially far from the cavity 
edge. The idea behind this set of experiments is to analyze both 
the orbital evolution of these new bodies, as well as their effects on 
the two coorbitals. Typical results are shown in Figure Q3] where 
each color curve corresponds to a different planet (see inlaid color 
code). Planets 1 and 2 (i.e. mi and rri2, respectively) are the origi- 
nal Saturn-mass bodies in coorbital motion, while planets 3 through 
8 are the new Earth-size masses. Since this type of simulation in- 
cluded smaller masses and a larger population, we employed a 
high-resolution description of the disk. 

We perceived no significant instability in the coorbitals, and 
their configuration remained unaffected by the new masses. These, 
however, suffered several cases of mutual scattering and close en- 
counters that caused the ejection of three of the masses. One of the 
close encounters (between planets 4 and 6) led to a temporary trap- 
ping of planet 4 in the cavity edge but finally led to the formation 
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Figure 16. Evolution of characteristics angles for coorbital solutions (top 
frame AA and middle frame Ara during the hydro-simulation. We reserve 
hand-left side to the system 1-2 and right side to the system 4-6. We plot 
in the lower frames the evolution of characteristics angles integrating the 
system without considering the interaction with the disk. At left frame two 
conditions were chosen, one at t=45() and the other at t= 11450 (marked as 
color circles in the middle frame), meanwhile at right bottom frame the final 
configuration of earth-like planets was integrated. 



of a new coorbital system with planet 8. The role of scattering be- 
tween planets as a formation mechanism for coorbitals has already 
been pointed out by Cresswell & Nelson (2006). The final outcome 
of this system consists of an inner coorbital system, a single Earth- 
mass planet trapped in an exterior mean-motion resonance and an 
additional coorbital pair in another resonance. Although this ap- 
pears a very complex multiple resonant configuration, we found no 
indication of long-term instability. 

Figure [T6] shows the dynamics of both pairs of coorbitals: 
planets 1 and 2 are displayed on the left-hand plots, while the coor- 
bital system composed of planets 4 and 8 are shown on the right. 
Top graphs present the temporal evolution of the resonant angle 
a = A A = A2 — Ai while the middle plots show Azu = zu2 — tn\. 
Once planet 4 is scattered into the cavity, the resonance relation 
between planets 1 and 2 is temporarily disrupted, resulting in a 
short-lived circulation of both angles. However, once the Earth- 
mass planet is sent back outside the cavity, the coorbital config- 
uration of the inner planets is reestablished, although around L5 in- 
stead of L4. We did some follow-up integrations without hydrody- 
namical interaction choosing two intermediate stages as initial con- 
ditions. The bottom left-hand plot shows the results choosing as ini- 
tial conditions the system configuration at t = 450 and t — 1 1450, 
respectively. The behavior of both critical angles show a small am- 
plitude libration around the equilibrium solutions. Between these 
two solutions the system experiments temporarily AL4, and AL5 
configurations. 

The dynamical evolution of the small-mass pair is more com- 
plex, and a stable configuration is only reached near the end of the 
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simulation, and corresponds to an AL5 type orbit with large ampli- 
tude of oscillation of Azu or about 150°. 



8 CONCLUSIONS 

In this paper we have analyzed the detectability and possible forma- 
tion mechanism of hypothetical massive planets in stable coorbital 
configurations. So far there are no known extrasolar planetary sys- 
tems containing coorbital bodies, which may imply that either these 
configurations are extremely difficult to form or to detect from RV 
surveys. 

We have studied the detectability of three different types of 
coorbital motion (QS, L4/L5 and AL4/AL5), trying to evaluate 
possible bias in detections and identify what kind of compatible 
configuration could be detected. The analysis of Keplerian contri- 
butions to radial velocities allowed us to predict the value for the 
signal of one single planet that could be confused with coorbital 
configurations. 

For even low values of errors in radial velocities measure- 
ments (e = 3 m/s) and observation time-spans covering four orbital 
periods (which include most of the presently detected exoplanets), 
coorbital configurations appear hard to identify: the results of the 
fitting process could easily confuse the RV data with that stemming 
from other configurations/systems (single planet or 2/1 MMR sys- 
tem). For large mass ratios, a correct identification of coorbital con- 
figuration is even more complicated and more easy to confuse with 
the other configurations. In all cases, the residuals of the different 
systems are comparable, even more so for large mass ratios (see 
Figg}. 

Observation data sets covering longer time-spans allow to de- 
tect mutual perturbations, but once again the best fits are not always 
associated to coorbital motion. We have found several cases where 
other resonant configurations (i.e. 3/1 commensurability) actually 
give smaller residuals and better results. 

Coorbital motion is not only sensitive to the data set, but also 
to the fitting procedure. Nested two-planet strategies may also con- 
fuse the real dynamics and yield results widely different from the 
nominal orbits; in this sense simultaneous two-planet fits seem 
more robust. 

Transit observations should help distinguish coorbitals planets 
from other solutions. 

Finally, we have found that stable coorbital systems with two 
massive planets may be formed from originally non-resonant orbits 
through their interaction with a inner cavity in the protoplanetary 
disk, as long as the surface density of the disk is sufficiently large. 
Both our N-body and hydro-simulations indicate a preference to- 
wards coorbitals with similar masses (i.e. mi ~ rri2). In all our 
simulations with large mass ratios, the smaller planet was either 
pushed inside the cavity or trapped in a mean-motion commensu- 
rability outside the density jump. 
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APPENDIX. ABOUT THE SCALE OF SIMULATIONS 

All our simulations were performed with a density jump centered 
at ro = 1. If we wish relate our results to anothor scale, for ex- 
ample ro = 0.01, we need to modify not only the spatial scale 
of the orbital system but also the surface density value. Here we 
present a simple recipe for the scale transformations in order to 
preserve the complete dynamics, including both the gravitational 
interactions and migration timescale. 

The equation of motion of a planet with mass m,: interacting 
gravitationally with other TV bodies of mass rrij is given by 

dt» - 2- r% 6lJ Uy) 

where are the unit vectors corresponding to the relative posi- 
tions. If we introduce a scale change in the coordinates defined by 
r = ar' and a temporal transformation defined by t = fit 1 , then 
the gravitational dynamics is invariant if = a 3 ^ 2 . 
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In our problem, however, the dynamical evolution of the plan- 
ets also includes their interaction with the gas disk, which is spec- 
ified by the gradient of the total differential torque Y. If we then 
desire to preserve the complete dynamics of the system under the 
spatial rescaling, then the torque must scale as: 

T = a _1 r'. (20) 

For a type I migration, the total torque is given by T(r) oc 
£(r)n 2 r 4 , where n is the mean motion. We also assume that the 
surface density is given by a power law expression of type: 

E(r) = E (-V (21) 



jo 

for a suitable exponent r\. The coefficient E denotes the surface 
density at r = ro. Then, if we apply the spatial and temporal trans- 
formations required to preserve the gravitational dynamics, we ob- 
tain 

r = aE f4Vn'V 4 = «r5fV'- (22) 



Therefore, for scale our simulation to a new suitable value of r () , it 
is necessary to do the transformation: 



r = ar (23) 
t = a 3/ Y (24) 
E = oT 2 T! (25) 



For this new radii, time and density value, the dynamics of system 
is invariant. 



